clear,clc,clf
N = [-3 -2 -1 0 1 3 3 2 5 6 7 6 9 11]; % 序号序列
X = [0 2 3 3 2 3 0 -1 -2 -3 -4 -5 1 2]; % 值序列
stem(N,X); % 绘制离散图
hold on;
plot(N,zeros(1,length(X)),'r') % 绘制横轴，zeros(1,n)为产生1行n列元素值为0的数组
set(gca,'box','on'); % 产生坐标轴设在方框上
xlabel('序列号');
ylabel('序列值');

%% 绘制原始信号及采样后的信号示例
clear,clc,clf
dt = 0.01;n = 0:90-1;
t = n*dt;
f = 10; % 原始信号频率10Hz
x = sin(3*pi*f*t+0.5); % 在计算机上的原始信号
dt = 0.1;
n = 0:10-1;
t1 = n*dt;
% 以10Hz的采样频率采样，取一样的采样长度
% 序号长度为原始信号序号长度的1/10
x1 = sin(3*pi*f*t1+0.5);
subplot(3,1,1);plot(t,x);
ylim([-1 1]);title('原始信号');
subplot(3,1,2);plot(t,x,t1,x1,'rp');
ylim([-1 1]);title('采样过程');
subplot(3,1,3);plot(t1,x1);
ylim([-1 1]);
xlabel('时间/s');
title('采样后信号');

%% 5.2 典型离散时间序列
% 5.2.1 单位冲激序列
% matlab中冲激序列可以用zeros函数实现
N = 5;
n = 1;
x =zeros(1,N); %产生N个点的单位冲激序列
x(n) = 1; % 第n个点的值为1

% 单位阶跃序列u(n)
% 在matlab中可以利用ones函数产生一个全1向量
% x = ones(1,N)
% 也可以利用逻辑关系表达式产生单位阶跃序列
% x = [(n-n0)>=0];
% 产生单位阶跃序列u(n),u(n-10)
n = 40;
x = ones(1,n);
xn = 0:n-1;
subplot(211);stem(xn,x);
grid on
axis([-1 51 0 1.1]);title('单位阶跃序列u(n)')
ylabel('u(n)');xlabel('n')
x = [zeros(1,10),1,ones(1,29)];
subplot(212);stem(xn,x);grid on
axis([-1 51 0 1.1]);title('单位阶跃序列u(n-10)')
ylabel('u(n-10)');xlabel('n')


%% 正弦序列
clear,clc
n = 0:59;
x = sin(pi/5*n);
stem(n,x);grid on
axis([0,40,-1.5,1.5]);title('正弦序列')
xlabel('n');ylabel('h(n');

%% 实指数序列
clear,clc,clf
n = 0:20;
a1 = 1.6;a2 = -1.6;
x1 = a1.^n;
x2 = a2.^n;
subplot(221);stem(n,x1,'filled');
grid on 
xlabel('n');
ylabel('h(n)');
title('x(n)=1.6^{n}')
subplot(222);stem(n,x2,'filled');
grid on 
xlabel('n');
ylabel('h(n)');
title('x(n)=-1.6^{n}')

%% 复指数序列
% x = exp((signma+j*omega)*n);

n = 0:50;
A = 3;a=-1/9;b=pi/5;
x = A*exp((a+1i*b)*n);
subplot(221);stem(n,real(x),'filled');grid on
axis([0 30 -4 4]);title('实部');xlabel('n');
subplot(222);stem(n,imag(x),'filled');grid on
axis([0 30 -4 4]);title('虚部');xlabel('n');
subplot(223);stem(n,abs(x),'filled');grid on
axis([0 30 -4 4]);title('模');xlabel('n');
subplot(224);stem(n,angle(x),'filled');grid on
axis([0 30 -4 4]);title('相角');xlabel('n');

%% 序列相加相乘
% x = x1+x2;
% x = x1.*x2;
% 序列值累加累积
n1 = [1 2 3 4];
sum(n1)
x1 = [2 0.5 0.9 2 2];
x = prod(x1)

%% 序列值累加与乘积
n1 = [1 2 3 4];
sum(n1)
x1=[2 0.5 0.9 2 2];
x=prod(x1)

%% 序列翻转与移位
% y=fliplr(x）
% ny=-fliplr(nx)

%% 连续信号的尺度变换
%矩形波的尺度变换
clear,clc
t=-4:0.001:4;
T=2;
f=rectpuls(t,T); 
ft=rectpuls(2*t,T);
subplot(211);plot(t,f);axis([-4,4,-0.5,1.5])
subplot(212);plot(t,ft);axis([-4,4,-0.5,1.5])
% 三角波尺度变换
t=-3:0.001:3;
ft=tripuls(t,4,0.5);
subplot(211);plot(t,ft)
ft=tripuls(3*t,4,0.5);
subplot(212);plot(t,ft)

%% 连续信号奇偶分解
% 对函数f(t)=cos(t+1)+t进行奇偶分解
syms t;
f=sym('cos(t+1)+t');
f1=subs(f,t,-t);
g=1/2*(f+f1);
h=1/2*(f-f1);
subplot(311);ezplot(f,[-8,8]);title('原信号')
subplot(312);ezplot(g,[-8,8]);title('偶分量')
subplot(313);ezplot(h,[-8,8]);title('奇分量')

%% 信号积分微分
% 微分运算用diff函数完成
% 调用格式 diff(function,'variable',n)
% function=需要进行求导运算的信号
% variable=求导运算的独立变量
% 连续运算的积分运算用int函数完成
% int(function,'variable',a,b)
% function=需要进行被积的信号
% variable=求积运算的独立变量
% a、b=积分上下限

% 微分&积分运算
syms t f2;
% heaviside
% this MATLAB function evaluates the Heaviside step 
% function (also known as the unit step function) at x.
f2=t*(2*heaviside(t)-heaviside(t-1))+heaviside(t-1);
t=-1:0.01:2;
subplot(121);ezplot(f2,t);grid on
ylabel('x(t)');title('原函数')
f=diff(f2,'t',1);
% f=int(f1,'t');
subplot(122);ezplot(f,t);grid on
ylabel('x(t)');title('微分后函数')

% 卷积运算
A=ones(1,3);
B=[1 8 8 5];
C=conv(A,B)
%% 5.4.1 随机信号函数
% rand(n) 返回n x n 的随机矩阵
% rand(m,n) 返回m x n 的随机矩阵
% rand(m,n,p,...) 产生随机数组
% rand(size(A)) 返回一个与A相同尺寸的随机矩阵

% 示例
tn=0:50;
N=length(tn);
x=rand(1,N);
subplot(121);plot(tn,x,'k');ylabel('x(t)');
subplot(122);stem(tn,x,"filled",'k');ylabel('x(n)');

%% 方波函数
% x=square(t)
% x=square(t,duty)
% 产生指定周期方波，duty为占空比
f=6000;nt=3;
N=15;T=1/f;dt=T/N;
n=0:nt*N-1;
tn=n*dt;
x=square(2*f*pi*tn,25)+1;
subplot(121);stairs(tn,x,'k');
axis([0 nt*T 1.1*min(x) 1.1*max(x)]);
ylabel('x(t)');
subplot(122);stem(tn,x,"filled",'k');
axis([0 nt*T 1.1*min(x) 1.1*max(x)]);
ylabel('x(n)');


%% 非周期方波
% rectpuls(t)
t=-3:0.001:3;
y=rectpuls(t);
subplot(121);plot(t,y);
axis([-2 2 -1 2]);grid on
xlabel('t');ylabel('w(t)');
y=2.5*rectpuls(t,2);
subplot(122);plot(t,y);
axis([-2 2 -1 3]);grid on
xlabel('t');ylabel('w(t)');

%% 三角波
Fs=10000;
t=0:1/Fs:1;
x1=sawtooth(2*pi*40*t,0);
x2=sawtooth(2*pi*40*t,1);
x3=sawtooth(2*pi*40*t,0.5);
subplot(311);plot(t,x1);axis([0 0.25 -1 1])
subplot(312);plot(t,x2);axis([0 0.25 -1 1])
subplot(313);plot(t,x3);axis([0 0.25 -1 1])

%% 非周期三角波
% tripuls
t=-2:0.001:2;
y=tripuls(t,4,0.5);
plot(t,y);grid on ;
axis([-3 3 -1 2]);
xlabel('t');
ylabel('y(t');

%% sinc函数
t=(1:12)';
x=randn(size(t));
ts=linspace(-10,10,500)';
y=sinc(ts(:,ones(size(t)))-t(:,ones(size(ts)))')*x;
plot(t,x,'o',ts,y);grid on 
xlabel('n')
ylabel('x(n)')

%% diric函数
t=-3*pi:pi/40:4*pi;
subplot(211);plot(t,sinc(t));grid on 
xlabel('t');ylabel('sinc(t)');title('Sinc');
subplot(212);plot(t,diric(t,5));grid on 
xlabel('t');ylabel('diric(t)');title('Diric');

%% 高斯调制正弦脉冲信号函数
% yi=gauspuls(t,fc,bw)
% yi=gauspuls(t,fc,bw,bwr)
% [yi,yq]=gauspuls()
% [yi,yq,ye]=gauspuls()
% tc=gauspuls('cutoff',fc,bw,bwr,tpe)

% 例5-26
tc=gauspuls('cutoff',60e3,0.6,[],-40);
t=-tc:1e-6:tc;
yi=gauspuls(t,60e3,0.6);
plot(t,yi);grid on
xlabel('t');ylabel('h(t')

%% 冲激串函数
T=0:1/1E3:1;
D=0:1/4:1;
Y=pulstran(T,D,'rectpuls',0.1);
subplot(211);plot(T,Y);axis([0 1 -0.1 1.1]);grid on
xlabel('t');ylabel('w(t)')
T=0:1/1E3:1;
D=0:1/3:1;
Y=pulstran(T,D,'tripuls',0.2,1);
subplot(212);plot(T,Y);axis([0 1 -0.1 1.1]);grid on
xlabel('t');ylabel('w(t)')


%% 线性调频函数
% chirp(t,f0,t,f1) 产生一个线性（频率随时间变化线性变化）信号，
% 其时间轴设置由数组t定义
% chirp(t,f0,t,f1，method) 指定改变扫频的方法
t=0:0.01:2;
y=chirp(t,0,1,120);
plot(t,y);axis([0 1 0 1]);grid on
ylabel('x(t');
xlabel('t');

%% 计算谱图与现行调频信号瞬时频率偏差
t=0:0.002:2;
y=chirp(t,0,1,150);
subplot(311);spectrogram(y,256,250,256,1E3,'yaxis');
xlabel('t=0:0.002:2');title('不同采样时间的条件下')
t=-2:0.002:2;
y=chirp(t,100,1,200,'quadratic');
subplot(323);spectrogram(y,128,120,128,1E3,'yaxis');
xlabel('t=-2:0.002:2');
t=-1:0.002:1;
fo=100;f1=400;
y=chirp(t,fo,1,f1,'q',[],'convex');
subplot(324);spectrogram(y,256,200,256,1000,'yaxis')
xlabel('t=-1:0.002:1');
t=0:0.002:1;
fo=100;f1=25;
y=chirp(t,fo,1,f1,'q',[],'concave');
subplot(325);spectrogram(y,hanning(256),128,256,1000,'yaxis')
xlabel('t=0:0.002:1');
t=0:0.002:10;
fo=10;f1=400;
y=chirp(t,fo,10,f1,'logarithmic');
subplot(326);spectrogram(y,256,128,256,1000,'yaxis')
xlabel('t=0:0.002:10');

% spectrogram函数是matlab中用于画频谱图的函数。
% 其中x是输入信号，
% Fs是采样率，
% window和noverlap是控制STFT(Synthesized Time-Frequency Transform)参数的变量，
% nfft是FFT计算的点数，
% fscale是控制颜色表缩放的变量。

% 其中，x与Fs是必填的参数，而window、noverlap、nfft和fscale则可以根据需要选择输入或使用默认参数。
% function [S,F,T] = spectrogram(x,window,noverlap,nfft,fs,validRange)
% S = spectrogram(x,window,noverlap,nfft,fs) returns the short-time Fourier
% transform of the input signal x. 
% 
% [S,F,T] = spectrogram(x,window,noverlap,nfft,fs) returns the STFT data S,
% frequency vector F and time vector T.

%% 高斯单脉冲函数
% y=gmonopuls(t,fc)
% 产生中心频率为fc的高斯单脉冲信号，t为采样间隔，默认fc=1000Hz

fc=2e9; %中心频率
fs=100e9; % 采样频率
tc=gmonopuls('cutoff',fc);  % 脉冲宽度
t=-2*tc:1/fs:2*tc;
y=gmonopuls(t,fc);
sg=1/(2*pi*fc);
ys=exp(1/2)*t/sg.*exp(-(t/sg).^2/2);
subplot(211);plot(t,y,t,ys,'.')
legend('gmonopuls','Definition')
D=((0:2)*7.5+2.5)*1e-9;
t=0:1/fs:150*tc;
yp=pulstran(t,D,'gmonopuls',fc);
subplot(212);plot(t,yp)

% pulstran 函数基于连续的或采样的原型脉冲生成脉冲序列。以下示例生成由高斯脉冲的
% 多次延迟插值之和组成的脉冲序列。该脉冲序列定义为具有 50 kHz 的采样率、10 ms 
% 的脉冲序列长度和 1 kHz 的脉冲重复率；D 在第 1 列中指定每个脉冲重复的延迟，
% 在第 2 列中指定每个重复的可选衰减。该脉冲序列是通过将 gauspuls 函数的名称
% 以及附加参数（用于指定带宽为 50% 的 10 kHz 高斯脉冲）传递给 pulstran 来构造的
%% 5.5 线性连续时间系统的时域分析
% 零状态与零输入的求解分析
syms y(t)
eqn=diff(y,2)+3*diff(y)+2*y==0; % 齐次解，求零输入响应
Dy=diff(y,t);
cond=[y(0)==1,Dy(0)==2];
yzi=dsolve(eqn,cond);
yzi=simplify(yzi)

%% 
ts=0;te=5;dt=0.01;
sys=tf(1,[1 2 200]);
t=ts:dt:te;
f=10*cos(2*pi*t);
y=lsim(sys,f,t);
plot(t,y);title('零状态响应');grid on
xlabel('t/s');ylabel('y(t)');

%% 冲激响应和阶跃响应分析
t=0:0.002:4;
sys=tf([1,32],[1,4,64]);
h=impulse(sys,t);
g=step(sys,t);
subplot(211);plot(t,h);title('冲激响应');grid on
xlabel('时间/s');ylabel('h(t)');
subplot(212);plot(t,g);title('阶跃响应');grid on
xlabel('时间/s');ylabel('g(t)');

%% 例5-34
clear,clc,clf
b=[-0.48 -0.25 -0.12 -0.06];a=[1 0.64 0.94 0.51 0.01];
sys=tf(b,a);
T=1000;
t=0:1/T:10;t1=-5:1/T:5;
f1=stepfun(t1,-1/T)-stepfun(t1,1/T);
f2=stepfun(t1,0);
f3=t;
f4=sin(t);
y1=lsim(sys,f1,t);y2=lsim(sys,f2,t);
y3=lsim(sys,f3,t);y4=lsim(sys,f4,t);
subplot(221);plot(t,y1);title('冲激激励下的零状态响应')
xlabel('t');ylabel('y1(t)');
grid on;axis([0 10 -1.2 1.2]);
subplot(222);plot(t,y2);title('阶跃激励下的零状态响应')
xlabel('t');ylabel('y2(t)');
grid on;axis([0 10 -1.2 1.2]);
subplot(223);plot(t,y3);title('斜坡激励下的零状态响应')
xlabel('t');ylabel('y3(t)');
grid on;axis([0 10 -5 0.5]);
subplot(224);plot(t,y4);title('正弦激励下的零状态响应')
xlabel('t');ylabel('y4(t)');
grid on;axis([0 10 -1.5 1.2]);

%% 卷积求解
dt=0.01;t=-1:dt:2.5;
f1=heaviside(t)-0.5*heaviside(t-2);
f2=2*exp(-3*t).*heaviside(t);
f=conv(f1,f2)*dt;n=length(f);tt=(0:n-1)*dt-2;
subplot(221);plot(t,f1);grid on 
axis([-1 2.5 -0.2 1.2]);title('f1(t)')
xlabel('t');ylabel('f1(t)');
subplot(222);plot(t,f2);grid on 
axis([-1 2.5 -0.2 1.2]);title('f2(t)')
xlabel('t');ylabel('f2(t)');
subplot(212);plot(tt,f);grid on 
xlabel('t');ylabel('f(t)');

%% 离散时间系统响应
a=[1 -0.35 1.5];
b=[1 1];
t=0:20;
x=(1/2).^t;
y=filter(b,a,x)
subplot(121);stem(t,x);title('离散序列');
xlabel('n');ylabel('h(n)');grid on
subplot(122);stem(t,y);title('响应序列');
xlabel('n');ylabel('h(n)');grid on
%% 离散时间系统的冲激响应和阶跃响应
k=0:10;
a=[1 64];
b=[1 3];
h=impz(b,a,k);
subplot(121);stem(k,h);title('单位冲激响应的近似值');
xlabel('n');ylabel('h(n)');grid on
hk=-(-1).^k+2*(-2).^k;
subplot(122);stem(k,h);title('单位冲激响应的理论值');
xlabel('n');ylabel('h(n)');grid on
%% 离散时间信号的卷积和运算
nx=-1:5;
nh=-2:10;
x=uDT(nx)-uDT(nx-4);
h=0.9.^nh.*(uDT(nh)-uDT(nh-9));
y=conv(x,h);
ny1=nx(1)+nh(1);
ny=ny1+(0:(length(nx)+length(nh)-2));
subplot(131);stem(nx,x,'fill');grid on
xlabel('n');ylabel('x(n)');title('x(n)')
axis([-4 16 0 2])
subplot(132);stem(nh,h,'fill');grid on
xlabel('n');ylabel('h(n)');title('h(n)')
axis([-4 16 0 2])
subplot(133);stem(ny,y,'fill');grid on
xlabel('n');ylabel('y(n)');title('y(n)=x(n)*h(n)')
axis([-4 16 0 4])


